{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Example: Static inverse free-boundary equilibrium calculations\n", "\n", "---\n", "\n", "This example notebook shows how to use FreeGSNKE to solve **static inverse** free-boundary Grad-Shafranov (GS) problems. \n", "\n", "In the **inverse** solve mode we seek to estimate the active poloidal field coil currents using user-defined constraints (e.g. on isoflux, x-point, o-point, and psi values) and plasma current density profiles for a desired equilibrium shape. \n", "\n", "Note that during this solve, currents are **not** found in any specified passive structures. \n", "\n", "Below, we illustrate how to use the solver for both diverted and limited plasma configurations in a **MAST-U-like tokamak** using stored pickle files containing the machine description. These machine description files partially come from the FreeGS repository and are not an exact replica of MAST-U." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The static free-boundary Grad-Shafranov problem\n", "\n", "Here we will outline the **static free-boundary** GS problem that is solved within both the **forward** and **inverse** solvers, though we encourage you to see [Pentland et al. (2024)](https://arxiv.org/abs/2407.12432) for more details. \n", "\n", "Using a cylindrical coordinate system $(R,\\phi,Z)$, the aim is to solve the GS equation:\n", "\n", "$$ \\Delta^* \\psi \\equiv \\left( \\frac{\\partial^2}{\\partial R^2} - \\frac{1}{R} \\frac{\\partial}{\\partial R} + \\frac{\\partial^2}{\\partial Z^2} \\right) \\psi = -\\mu_0 R J_{\\phi}(\\psi, R, Z), \\qquad (R,Z) \\in \\Omega, $$\n", "\n", "for the poloidal flux $\\psi(R,Z)$ (which here has units Weber/$2\\pi$) in the rectangular computational domain $\\Omega$. The flux has contributions from both the plasma and the coils (metals) such that $\\psi = \\psi_p + \\psi_c$. This flux defines the toroidal current density $J_{\\phi} = J_p(\\psi,R,Z) + J_c(R,Z)$, also containing a contribution from both the plasma and coils, respectively. We have the plasma current density (only valid in the core plasma region $\\Omega_p$):\n", "\n", "$$ J_p(\\psi,R,Z) = R \\frac{\\mathrm{d}p}{\\mathrm{d}\\psi} +\\frac{1}{\\mu_0 R} F \\frac{\\mathrm{d} F}{\\mathrm{d} \\psi}, \\qquad (R,Z) \\in \\Omega_p, $$\n", "\n", "where $p(\\psi)$ is the plasma pressure profile and $F(\\psi)$ is the toroidal magnetic field profile. The current density generated by $N$ active coils and passive structures is given by:\n", "\n", "$$ J_c(R,Z) = \\sum_{j=1}^{N} \\frac{I^c_j(R,Z)}{A_j^c}, \\qquad (R,Z) \\in \\Omega, \\quad \\text{where} \\quad I_j^c(R,Z) = \n", "\\begin{cases} \n", " I_j^c & \\text{if } (R,Z) \\in \\Omega_j^c, \\\\ \n", " 0 & \\text{elsewhere}.\n", "\\end{cases}.$$\n", "\n", "This makes use of the current $I^c_j$ in each metal and its cross-sectional area $A^c_j$ (the domain of each metal is denoted by $\\Omega_j^c$).\n", "\n", "To complete the problem, we need the integral (Dirichlet) free-boundary condition:\n", "\n", "$$ \\psi(R,Z) = \\int_{\\Omega} G(R,Z;R',Z') J_{\\phi}(\\psi, R',Z') \\ \\mathrm{d}R' \\mathrm{d}Z', \\qquad (R,Z) \\in \\partial \\Omega, $$\n", "\n", "where $G$ is the (known) Green's function for the elliptic operator above.\n", "\n", "We'll now go through the steps required to solve the **inverse** problem in FreeGSNKE. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Create the machine object\n", "\n", "First, we build the machine object from previously created pickle files in the \"machine_configs/MAST-U\" directory. \n", "\n", "FreeGSNKE requires the following paths in order to build the machine:\n", "- `active_coils_path`\n", "- `passive_coils_path`\n", "- `limiter_path`\n", "- `wall_path`\n", "- `magnetic_probe_path` (not required here)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# build machine\n", "from freegsnke import build_machine\n", "tokamak = build_machine.tokamak(\n", " active_coils_path=f\"../machine_configs/MAST-U/MAST-U_like_active_coils.pickle\",\n", " passive_coils_path=f\"../machine_configs/MAST-U/MAST-U_like_passive_coils.pickle\",\n", " limiter_path=f\"../machine_configs/MAST-U/MAST-U_like_limiter.pickle\",\n", " wall_path=f\"../machine_configs/MAST-U/MAST-U_like_wall.pickle\",\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot the machine\n", "import matplotlib.pyplot as plt\n", "\n", "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)\n", "plt.tight_layout()\n", "\n", "tokamak.plot(axis=ax1, show=False)\n", "ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "\n", "ax1.grid(alpha=0.5)\n", "ax1.set_aspect('equal')\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", "ax1.set_ylabel(r'Height, $Z$ [m]')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Instantiate an equilibrium\n", "\n", "We are now ready to build a plasma equilibrium object for our tokamak. This is done using the `freegs4e.Equilibrium` class, which implicitly defines the rectangular domain of the solver as well as the grid resolution.\n", "\n", "`Equilibrium` has sensible defaults, but it is recommended to define the radial and vertical domain of the grid using the `Rmin`, `Rmax`, `Zmin` and `Zmax` parameters, as well as the grid resolution in the radial and vertical directions with the `nx` and `ny` parameters. The grid will be initialised using fourth-order finite differences. Note that the computational grid **must** encompass the domain enclosed by the limiter object (as this is where the plasma will be confined to). \n", "\n", "A tokamak object should be supplied to the `tokamak` parameter to assign the desired machine to the equilibrium.\n", "\n", "If available, an initial guess for the plasma flux $\\psi_p$ (dimensions `nx` x `ny`) can be provided via the `psi` parameter (commented out in the following code). If not, the default initialisation will be used. \n", "\n", "The `eq` object will store a lot of important information and derived quantites once the equilibrium has been calculated (see future notebook on this). " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from freegsnke import equilibrium_update\n", "\n", "eq = equilibrium_update.Equilibrium(\n", " tokamak=tokamak, # provide tokamak object\n", " Rmin=0.1, Rmax=2.0, # radial range\n", " Zmin=-2.2, Zmax=2.2, # vertical range\n", " nx=65, # number of grid points in the radial direction (needs to be of the form (2**n + 1) with n being an integer)\n", " ny=129, # number of grid points in the vertical direction (needs to be of the form (2**n + 1) with n being an integer)\n", " # psi=plasma_psi\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Instantiate a profile object\n", "\n", "We can now instatiate a profile object that contains the chosen parameterisation of the toroidal plasma current density $J_p$ (i.e. on right hand side of the GS equation). We can then set the paramters for the chosen current density profiles. \n", "\n", "A number of commonly used profile parameterisations exist in FreeGSNKE, including:\n", "- `ConstrainPaxisIp`\n", "- `ConstrainBetapIp`\n", "- `Fiesta_Topeol`\n", "- `Lao85`\n", "- `TensionSpline`\n", "\n", "In this notebook, we will make use of the `ConstrainPaxisIp` (and `ConstrainBetapIp`) profiles (see [Jeon (2015)](https://link.springer.com/article/10.3938/jkps.67.843)). Others will be utilised in later notebooks. If there is a profile parameterisation you require that does not exist, please do create an issue. \n", "\n", "Both `ConstrainPaxisIp` and `ConstrainBetapIp` are parameterised as follows:\n", " $$J_{p}(\\psi, R, Z) = \\lambda\\big[ \\beta_{0} \\frac{R}{R_{0}} \\left( 1-\\tilde{\\psi}^{\\alpha_m} \\right)^{\\alpha_n} + (1-\\beta_{0}) \\frac{R_0}{R} \\left( 1-\\tilde{\\psi}^{\\alpha_m} \\right)^{\\alpha_n} \\big] \\quad (R,Z) \\in \\Omega_p,$$\n", "\n", "where the first term is the pressure profile and the second is the toroidal current profile. Here, $\\tilde{\\psi}$ denotes the normalised flux:\n", "$$ \\tilde{\\psi} = \\frac{\\psi - \\psi_a}{\\psi_b - \\psi_a}, $$\n", "where $\\psi_a$ and $\\psi_b$ are the values of the flux on the magnetic axis and plasma boundary, respectively. \n", "\n", "The parameters required to define this particular profile are:\n", "- `Ip` (total plasma current).\n", "- `fvac` ($rB_{tor}$, vacuum toroidal field strength).\n", "- `alpha_m`>0, and `alpha_n`>0 (that define the shape/peakedness of the profiles).\n", "- If `ConstrainPaxisIp` is used, then `paxis` (pressure on the magnetic axis) is required.\n", "- If `ConstrainBetapIp` is used, then `betap` (proxy of the poloidal beta) is required.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The values of $\\lambda$ and $\\beta_0$ are found using the above parameters as constraints ($R_0$ is a fixed scaling constant) in the following:\n", "\n", "- For `ConstrainPaxisIp`, we can re-arrange the following equations to solve for the unknowns:\n", "\n", "$$ p_{\\text{axis}} = \\lambda \\beta_{0} \\frac{R}{R_{0}} \\int^{\\psi_b}_{\\psi_a} \\left( 1-\\tilde{\\psi}^{\\alpha_m} \\right)^{\\alpha_n} \\mathrm{d} \\tilde{\\psi} $$\n", "\n", "and\n", "\n", "$$ I_p = \\int^{Z_{\\text{max}}}_{Z_{\\text{min}}} \\int^{R_{\\text{max}}}_{R_{\\text{min}}} J_p(\\psi, R, Z) \\ \\mathrm{d}R \\mathrm{d}Z. $$\n", "\n", "\n", "- For `ConstrainBetapIp`, we can instead re-arrange and solve the following:\n", "\n", "$$ \\beta_p = \\frac{8 \\pi}{\\mu_0 I_p^2} \\int^{Z_{\\text{max}}}_{Z_{\\text{min}}} \\int^{R_{\\text{max}}}_{R_{\\text{min}}} p(\\psi) \\ \\mathrm{d}R \\mathrm{d}Z. $$\n", "\n", "and\n", "\n", "$$I_p = \\int^{Z_{\\text{max}}}_{Z_{\\text{min}}} \\int^{R_{\\text{max}}}_{R_{\\text{min}}} J_p(\\psi, R, Z) \\ \\mathrm{d}R \\mathrm{d}Z. $$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In what follows, we use `ConstrainPaxisIp`. Note that the equilibrium (`eq`) object is passed to the profile to inform calculations relating to the machine description." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# initialise the profiles\n", "from freegsnke.jtor_update import ConstrainPaxisIp\n", "profiles = ConstrainPaxisIp(\n", " eq=eq, # equilibrium object\n", " paxis=8e3, # profile object\n", " Ip=6e5, # plasma current\n", " fvac=0.5, # fvac = rB_{tor}\n", " alpha_m=1.8, # profile function parameter\n", " alpha_n=1.2 # profile function parameter\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Load the static nonlinear solver\n", "\n", "We can now load FreeGSNKE's Grad-Shafranov static solver. The equilibrium is used to inform the solver of the computational domain and of the tokamak properties. The solver below can be used for **both** inverse and forward solves. \n", "\n", "Note: It's not necessary to instantiate a new solver when aiming to use it on new or different equilibria, as long as the integration domain, mesh grid, and tokamak are consistent across solves. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from freegsnke import GSstaticsolver\n", "GSStaticSolver = GSstaticsolver.NKGSsolver(eq) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Constraints\n", "\n", "Recall that in the **inverse** solve mode we seek to **estimate the active poloidal field coil currents** using user-defined **constraints** (e.g. on isoflux, null points, and psi values) and plasma current density profiles for a desired equilibrium shape. \n", "\n", "FreeGSNKE uses a `constrain` object, which accepts constraints on the location of:\n", "- any null points, which are either X-points or O-points (`null_points`).\n", "- any pairs of points that lie on the same flux surface (`isoflux_set`). The flux at these locations will be constrained to be the same. Multiple sets can be defined. \n", "- any known flux points (`psi_vals`), i.e. one can specify $\\psi$ at any location $(R,Z)$, possibly an entire flux map (not used here). \n", "\n", "At least one constraint (preferably many more) is required to carry out an inverse solve.\n", "\n", "Here, we specify two X-point locations (we want a double null plasma) and a number of isoflux locations. The isofluxes here will define the core plasma shape and part of the divertor legs. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np \n", "\n", "Rx = 0.6 # X-point radius\n", "Zx = 1.1 # X-point height\n", "Ra = .85\n", "Rout = 1.4 # outboard midplane radius\n", "Rin = 0.34 # inboard midplane radius\n", "\n", "# set desired null_points locations\n", "# this can include X-point and O-point locations\n", "null_points = [[Rx, Rx], [Zx, -Zx]]\n", "\n", "# set desired isoflux constraints with format \n", "# isoflux_set = [isoflux_0, isoflux_1 ... ] \n", "# with each isoflux_i = [R_coords, Z_coords]\n", "isoflux_set = np.array([[[Rx, Rx, Rin, Rout, 1.3, 1.3, .8,.8], [Zx, -Zx, 0.,0., 2.1, -2.1,1.62,-1.62]]])\n", " \n", "# instantiate the freegsnke constrain object\n", "from freegsnke.inverse import Inverse_optimizer\n", "constrain = Inverse_optimizer(null_points=null_points,\n", " isoflux_set=isoflux_set)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Given that there may be more or less constraints than unknown parameters (i.e. coil currents), the inverse problem may be over- or under-constrained. This means that there may be zero, one, or many solutions to the problem. \n", "\n", "A number of quadratic (i.e.Tikhonov) regularization parameters will be used to combat this. Larger values will encourage lower absolute coil current values. It is sometimes useful to experiment with different values to explore whether the converged solution departs from the desired constraints. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The linear system of constraints\n", "\n", "During an inverse solve, a minmisation problem involving the responses, changes in coil currents, and constraints, is repeatedly solved:\n", "\n", "$$ \\min_{x} \\| A\\vec{x} - \\vec{b}\\|^2 + \\| \\gamma \\vec{x} \\|^2, $$\n", "\n", "where\n", "- $A$ is the fixed response matrix (that determines how a change in coil currents $x$, affects constraint values). \n", "- $\\vec{x} = \\Delta \\vec{I}^c$ is the step change in active coil currents required to match the constraints. \n", "- $\\vec{b}$ is the vector of constraint values being enforced. \n", "- $\\gamma > 0$ is the regularisation parameter/vector.\n", "\n", "We solve for $\\vec{x}$ using a gradient-based optimiser. \n", "\n", "[Song et al. (2024)](https://www.mdpi.com/2571-6182/7/4/45) provide a nice overview of the inverse problem.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Fixed coil currents\n", "\n", "It is also possible to set values for the current in specific active poloidal field coils. The inverse solver will not allow the current in these coils to vary if the `control` parameter is set to `False` (i.e. it will be excluded during the optimisation). \n", "\n", "Note: any passive structures in the tokamak automatically have their control parameter set to False and are therefore not included in an inverse solve. \n", "\n", "As an example, we will fix the `Solenoid` current and seek a solution in which this value is fixed, rather than estimated by the inverse solve." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "eq.tokamak.set_coil_current('Solenoid', 5000)\n", "eq.tokamak['Solenoid'].control = False # ensures the current in the Solenoid is fixed" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The inverse solve\n", "\n", "The following cell will execute the solve. Since a `constrain` object is provided, this is interpreted as a call to the inverse solver, if `constrain=None`, then the forward solver will be called (see next notebook). The `target_relative_tolerance` is the maximum relative error on the plasma flux function allowed for convergence and `target_relative_psit_update` ensures that the relative update to the plasma flux (caused by the update in the control currents) is lower than this target value. Both are required to be met for the inverse problem to be considered successfully solved.\n", "\n", "The `verbose=True` option will provide an indication of the progression of the solve. \n", "\n", "The solver also needs stabilisation for any up-down symmetric coils used for vertical control of the plasma (in this case the `P6` coil). To do this, one can set the `l2_reg` parameter: the Tikonov regularisation used by the optmiser of the coil currents. Here, `P6` is the last coil in the list of those available for control, and so we use a larger regularisation for it specifically. This is the length of the coils being controlled (so it excludes `Solenoid`). \n", "\n", "The solver steps are (roughly):\n", "1. Solve the linear system to find initial coil currents that approximately satisfy the constraints for the initial equilibrium. \n", "2. While tolerance is not met:\n", " - use the coil currents to solve the forward GS problem (with NK iterations).\n", " - solve the linear system to update coil currents to satisfy constraints for current equilibrium (from forward solve). \n", " - check tolerances. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "GSStaticSolver.solve(eq=eq, \n", " profiles=profiles, \n", " constrain=constrain, \n", " target_relative_tolerance=1e-6,\n", " target_relative_psit_update=1e-3,\n", " verbose=True, # print output\n", " l2_reg=np.array([1e-12]*10+[1e-6]), \n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following plots show how to display:\n", "1. the tokamak with:\n", " - active coil filaments (rectangles with blue interior)\n", " - passive structures (blue circles if defined as filaments or thin black outline/grey interior if defined as parallelograms)\n", " - limiter/wall (solid black)\n", "2. the tokamak + the equilibrium with:\n", " - separatrix/last closed flux surface (solid red lines)\n", " - poloidal flux (yellow/green/blue contours, colours indicates magnitude)\n", " - X-points (red circles)\n", " - O-points (green crosses)\n", "3. the tokamak + the equilibrum + the contraints with:\n", " - null-point constraints (purple Y-shaped crosses)\n", " - isoflux contour constraints (usual crosses)\n", "\n", "Note: Setting 'show=True' can toggle the legend on/off." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig1, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12, 8), dpi=80)\n", "\n", "ax1.grid(zorder=0, alpha=0.75)\n", "ax1.set_aspect('equal')\n", "eq.tokamak.plot(axis=ax1,show=False) # plots the active coils and passive structures\n", "ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) # plots the limiter\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "\n", "ax2.grid(zorder=0, alpha=0.75)\n", "ax2.set_aspect('equal')\n", "eq.tokamak.plot(axis=ax2,show=False) # plots the active coils and passive structures\n", "ax2.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) # plots the limiter\n", "eq.plot(axis=ax2,show=False) # plots the equilibrium\n", "ax2.set_xlim(0.1, 2.15)\n", "ax2.set_ylim(-2.25, 2.25)\n", "\n", "\n", "ax3.grid(zorder=0, alpha=0.75)\n", "ax3.set_aspect('equal')\n", "eq.tokamak.plot(axis=ax3,show=False) # plots the active coils and passive structures\n", "ax3.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) # plots the limiter\n", "eq.plot(axis=ax3,show=False) # plots the equilibrium\n", "constrain.plot(axis=ax3, show=True) # plots the contraints\n", "ax3.set_xlim(0.1, 2.15)\n", "ax3.set_ylim(-2.25, 2.25)\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A solve call will modify the equilibrium object in place. That means that certain quantities within the object will be updated as a result of the solve. \n", "\n", "Various different quantities and functions can be accessed via the 'eq' and 'profile' objects. For example:\n", "- the total flux can be accessed with `eq.psi()`.\n", "- the plasma flux with `eq.plasma_psi`.\n", "- the active coil + passive structure flux with `eq.tokamak.psi(eq.R, eq.Z)`.\n", "- (Total flux = plasma flux + coil flux)\n", "\n", "Explore `eq.` to see more (also `profiles.`, e.g. the plasma current distribution over the domain can be found with `profiles.jtor`). Also see the `example3` notebook for how to access many different equilibrium-derived quantities of interest. \n", "\n", "The set of optimised coil currents can be retrieved using `eq.tokamak.getCurrents()` having been assigned to the equilibrium object during the inverse solve.\n", "\n", "The following lines will save the calculated currents to a pickle file (we will use these in future notebooks). " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "inverse_current_values = eq.tokamak.getCurrents()\n", "\n", "# save coil currents to file\n", "import pickle\n", "with open('simple_diverted_currents_PaxisIp.pk', 'wb') as f:\n", " pickle.dump(obj=inverse_current_values, file=f)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Second inverse solve: limited plasma example\n", "\n", "Below we carry out an inverse solve seeking coil current values for a limited plasma configuration (rather than a diverted one). \n", "\n", "In a limiter configuration the plasma \"touches\" the limiter of the tokamak and is confined by the solid structures of the vessel. The last closed flux surface (LCFS) is the closed contour that is furthest from the magnetic axis that just barely touches (and is tangent to) the limiter.\n", "\n", "We instantiate a new equilibrium but use the same profile (illustrating how to vary the total plasma current and pressure on axis within it) and solver objects." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "# a new eq object resets the coil currents and equilibrium\n", "eq = equilibrium_update.Equilibrium(\n", " tokamak=tokamak, # provide tokamak object\n", " Rmin=0.1, Rmax=2.0, # radial range\n", " Zmin=-2.2, Zmax=2.2, # vertical range\n", " nx=65, # number of grid points in the radial direction (needs to be of the form (2**n + 1) with n being an integer)\n", " ny=129, # number of grid points in the vertical direction (needs to be of the form (2**n + 1) with n being an integer)\n", " # psi=plasma_psi\n", ")\n", "\n", "# modify the total plasma current\n", "profiles.Ip = 4e5\n", "\n", "# modify the pressure on the magnetic axis\n", "profiles.paxis = 6e3\n", "\n", "# first we specify some alternative constraints\n", "Rout = 1.4 # outboard midplane radius\n", "Rin = 0.24 # inboard midplane radius\n", "\n", "# locations of X-points\n", "Rx = 0.45\n", "Zx = 1.18\n", "null_points = [[Rx, Rx], [Zx, -Zx]]\n", "\n", "# isoflux constraints\n", "isoflux_set = [[Rx, Rx, Rout, Rin, Rin, Rin, Rin, Rin, .75, .75, .85, .85, 1.3, 1.3],\n", " [Zx, -Zx, 0, 0, .1, -.1, .2, -.2, 1.6, -1.6, 1.7, -1.7, 2.1, -2.1]]\n", "\n", "# let's assume we're also seeking an equilibrium with no solenoid current\n", "eq.tokamak.set_coil_current('Solenoid', 0)\n", "eq.tokamak['Solenoid'].control = False # fixes the current\n", "\n", "# pass the magnetic constraints to a new constrain object\n", "constrain = Inverse_optimizer(null_points=null_points,\n", " isoflux_set=isoflux_set)\n", "\n", "# carry out the solve\n", "GSStaticSolver.solve(eq=eq, \n", " profiles=profiles, \n", " constrain=constrain, \n", " target_relative_tolerance=1e-6,\n", " target_relative_psit_update=1e-3,\n", " verbose=True, # print output\n", " l2_reg=np.array([1e-9]*10+[1e-5]), \n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# save the currents for later use\n", "inverse_current_values = eq.tokamak.getCurrents()\n", "\n", "# save coil currents to file\n", "with open('simple_limited_currents_PaxisIp.pk', 'wb') as f:\n", " pickle.dump(obj=inverse_current_values, file=f)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot the resulting equilbria \n", "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)\n", "ax1.grid(True, which='both')\n", "eq.plot(axis=ax1, show=False)\n", "eq.tokamak.plot(axis=ax1, show=False)\n", "constrain.plot(axis=ax1,show=True)\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "plt.tight_layout()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "freegsnke4e", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.16" } }, "nbformat": 4, "nbformat_minor": 4 }